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ABSTRACT 



As efforts continue to effectively utilize the microwave frequencies 
a need is apparent for a traveling wave type device capable of handling 
greater power at higher frequencies. One possible solution is the 
plasma microwave amplifier tube. 

This paper discusses an analysis of the plasma - helix coupling 
problem which may hopefully lead to development of such a device 
whereby the R. F. energy is coupled directly to the plasma with a 
helix, rather than coupling to an electron beam as has been the case 
heretofore. A theoretical analysis is presented, and the computer 
solution to the analytical equations evolved and the resulting W-/3 
plots are included. 

This writer wishes to express his appreciation for the assistance 
and inspiration provided by Dr. Glen A. Gray of the U. S. Naval 
Postgraduate School in carrying out this investigation. 



ii 



TABLE OF CONTENTS 



Page 

1 . Background 1 

2. Statement of the problem 11 

3. Analytical solution 12 

4. Discussion 19 

5. Conclusions 3 4 

6. Proposed Experimental Work 35 

7 . Bibliography 39 

8. Appendix I - computer program 41 



iii 



LIST OF ILLUSTRATIONS 



Page 



Fig. 


(1) 


Experimental Plasma Microwave Amplifier tube 


2 


Fig. 


(2) 


Experimental tube of Boyd, Field, and Gould (9) 


5 


Fig. 


(3) 


Experimental tube of Allen and Kino (8) 


6 


Fig. 


(4) 


Experimental tube of Paik (12) 


8 


Fig. 


(5) 


Sketch of Backward wave coupling 


10 


Fig. 


(6) 


Sketch of helix - plasma relationship showing con- 
ventions used 


13 


Fig. 


(7) 


co - # plot 

k a cot 0 = 3,0, b/a = 0,8, k a cot 0 /c = 5.0 
p c T 


23 


Fig. 


(8) 


Plot showing relationships of eq. (1) resulting in 
multiple modes of operation 


24 


Fig. 


(9) 


co- j3 plot 

k a cot 0 = 3.0, b/a = 0.8, k a cot 0 = 3.0 
P c 


25 


Fig. 


(10) 


CO - P plot 

k a cot 0=3.0, b/a = 0.8, k a cot 0=2.0 
P c 


26 


Fig. 


(ID 


co - /3plot 

k a cot 0= 3.0, b/a = 0.8, k a cot 0= 1.0 
P c 


27 


Fig. 


(12) 


co - 0 plot 

k a cot 0=3.0, b/a = 0.8, k a cot 0=0.5 
p c 


28 


Fig. 


(13) 


co - jSplot 

k a cot 0=3.0, b/a = 0.8, k a cot 0=0.1 
P c 


29 


Fig. 


(14) 


co - £ plot 

k a cot 0=3.0, b/a = 0.5, k a cot 0 /c = 5.0 
P c 


30 


Fig. 


(15) 


co - plot 

k a cot 0=3.0, b/a = 0.5, k a cot 0= 3.0 
P c 


31 



iv 



Fig. (16) 


oo - j8 plot 

k a cot i p = 3.0, b/a = 0 . 5 , k a cot ip = 2.0 
P c 


Page 

32 


Fig. (17) 


oo - jSplot 

k a cot ip = 3.0, b/a = 0 . 8 , k =0.0 
P c 

k a cot ip = 3.0, b/a = 0 . 8 , k -* 00 
p c 


33 


Fig. (18) 


Laboratory equipment 


37 


Fig. (19) 


Diagram of proposed experimental setup 


38 



V 



1. Background. 



The purpose of this investigation was to study the problem of 
coupling electromagnetic energy between a slow wave structure and a 
plasma column. 

The underlying reason for interest in this problem is the design 
of a plasma microwave amplifier tube. A plasma microwave amplifier 
tube, where the plasma column acts as the slow wave structure, inter- 
acting with the electron beam, may well provide some advantages over 
present day types of traveling wave tubes. At the present time, con- 
ventional traveling wave tubes are limited in the power they can handle 
at high microwave frequencies. As frequency increases, it becomes 
difficult to construct a satisfactory slow wave structure which will 
dissipate any appreciable heat. Even at the lower microwave frequencies, 
heat dissipation in the slow wave structure may be a problem for high 
power applications. The plasma as a slow wave structure offers 
definite possibilities in that heat dissipation in the structure should no 
longer be a problem. There should be no inherent limitation insofar as 
mechanical dimensions are concerned, as the frequency is increased. 
Further, the plasma microwave amplifier device appears to offer 
promise of high gain and high power in a relatively small package. 

The device proposed (Fig. (1) ) is in many respects similar to a 
conventional traveling wave tube. The tube consists of a conventional 
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Fig. 



electron gun assembly, a d.c. heated filament consisting of a .001" 
tungsten ribbon approximately 1/2" X 3/4" at either end of the interaction 
region, and two helix slow wave structures in the interaction region. 

The concept is to generate a plasma by heating the two filaments 
to 2300°C in a cesium gas environment at a very low pressure, on the 

-6 -7 

order of 10 to 10 mm. Hg. These hot filaments will then cause 
contact ionization of the cesium gas thereby generating a plasma. 

R. F. energy is then coupled from the helix slow wave structure to the 
plasma, exciting the plasma. In turn, the electron beam on passing 
thru the plasma interacts with the wave excited on the plasma, giving 
a traveling wave type interaction. 

A helix slow wave structure was chosen for several reasons. First, 
it is rather easily fabricated. Second, it has been used extensively 
in tube work, and experience is available regarding its characteristics,. 
Finally, a theoretical analysis of the helix is available. Beck (1), and 
H utter (2). 

The cesium plasma was chosen for several reasons. First, the 
plasma frequency may be easily varied by control of the envelope 
temperature, and thereby the gas pressure. Second, it provides a 
relatively homogenous plasma, and may be quite highly ionized. Third, 
it provides a quiet plasma, that is, it is relatively free from spurious 
oscillations and noise. Finally, it may be readily ionized by contact 
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ionization. Since the ionization voltage of the cesium is below the 
work function of the heated tungsten filament, the filament will emit 
electrons and at the same time ionize the cesium gas, generating a 
plasma which is relatively homogenous. Knechtli and Wada (7) have 
done considerable work on this type of plasma generation. Kino (8) 
has also used a similar scheme. In his tube. Fig. (3), a spiral heater 
was used, and due to the fields set up across the spiral, it was 
necessary to make all measurements during the time when power was 
not applied to the heater. By use of a flat ribbon, it is hoped to be 
able to reduce the field of the heater sufficiently so that measurements 
may be taken continously while applying a d.c. voltage to the heater. 

Work on plasma microwave amplifier tubes has been carried on 
by several authors , (8), (9), (12), (13), (14), and (15). The major 
difference between the operation of the tube proposed here and the 
work of some previous authors is in the manner of coupling the energy 
into the device. Boyd, Field, and Gould, (9) studied plasma oscillations 
and growing plasma waves using a mercury arc discharge for a plasma, 
coupling the input R. F. signal to the electron beam, and thence to the 
plasma. Fig. (2) is a sketch of the device they used. Allen and Kino (8) 
studied the interaction of an electron beam with a fully ionized plasma 
using a cesium plasma and coupling tie input R. F. signal to the electron 
beam, and thence to the plasma. Fig. (3) is a sketch of the device they 
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used. Paik (12) used a coupled mode approach to study the coupling of 
modes between a slow wave mode and a helix. The experimental work 
he reported on was done with a device as shown in Fig. (4) which used 
a mercury plasma. He excited only one mode by careful coupling, 
using a selective mode coupler. This device did not include an electron 
beam, and his theoretical and experimental work was limited to the case 
of no magnetic field. 

B.M. Bulgakov, et al. (13) studied slow waves in a helix - plasma 
column with an arbitrary magnetic field. This did not include any ex- 
perimental work, and the final dispersion plots they developed showed 
discontinuities in the curves, which one would not expect to see. 

This has caused some workers in the field to question their results. 

E. V. Bogdanov, etal. (14) studied very nearly the same problem 
as Boyd, Field, and Gould (9). The device they used for experimental 
work, although physically different, used an arc discharge type of 
plasma generator and coupled the input R. F. energy to an electron 
beam, and the output R. F. energy from the electron beam. The signifi- 
cant thing about this work is they were able to obtain gain from 25 cm. 
up to 3 cm. This was possible because they used an arc discharge to 
generate the plasma, and by increasing the arc current could generate 
a dense plasma, which was then contained by a magnetic field. 

Work on an experimental X band tube has been carried out by 
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Microwave Associates , Burlington, Massachusetts. 1 This is a mercury 
device, again coupling the signal via an electron beam. Eicon Labor- 
atory, Cambridge, Massachusetts is working on a device using the fast 
wave, designed to utilize the negative conductance of the plasma. All 
other devices proposed to date have used the slow wave mechanism. 
General Electric Ltd. , Wembly England, has reported* a plasma X band 
unit which develops 60 db of gain in a six cm. tube. Philips Research 
Laboratories, Eindhoven, Netherlands has reported* a helix coupled 
tube giving 19 db gain at four Gc. 

In all of the preceding devices coupling of the input signal to the 
device was accomplished by coupling as in a conventional traveling 
wave tube, to an electron beam with a helix, and then using the wave 
propagated on the beam to couple to the plasma. As a result of this, 
in the experimental work performed to date it has not been possible to 
couple directly to a backward wave plasma mode. It is anticipated that 
with the experimental device proposed. Fig. (1), it will be possible to 
investigate backward wave modes by direct coupling, as shown in 
Fig. (5). 



*M. Meisels, Plasma Amplifiers Move Closer to Hardware Stage, 
Microwaves, pp 4-6, Jan. 1963 
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If energy is coupled into and out of the plasma in this fashion, then _if_ 
interaction with the electron beam is observed, it must necessarily be 
interaction with a backward wave. Of course, care must be exercised 
to ascertain whether this interaction occurs in the helix or in the drift 
region. 

In considering the specific problem of coupling electromagnetic 
energy to a plasma column, it is perhaps worthwhile to review the work 
done in studying propagation on an electron beam. Hahn (10) and 
Ramo (11) considered a drifting, ion-neutralized cylindrical electron 
beam in an infinite axial magnetic field and found that in general there 
are two space charge waves, the propagation equation being. 

to 2 - (cc - 0uo) 2 = 0 (1.1) 

P 

for the one dimensional case, where 
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( 1 . 2 ) 




and all quantities are assumed to have e^ ^ 1 ^ z dependence. Thus, 
the solutions are; 



fi = 



Uo 



(1.3) 



If Uq is zero, this one dimensional ion- neutralized electron beam 
does not propagate space charge disturbances. Trivelpiece (6) has 
shown however, that a finite u 0 is not essential to the propagation of 
space charge disturbances when the beam has a finite transverse geo- 
metry. 

Rigrod and Lewis (3) considered the problem of an electron beam 
under Brillouin flow conditions and found two types of space charge 
wave propagation, one involving a perturbation of the average charge 
density, and the other surface rippling. Brewer (4) considered the more 
general case of an arbitrary magnetic field. 



2. Statement of the problem. 

In order to study the coupling of R.F. energy between a helix slow 
wave structure and a plasma column, it is desirable to determine the 
CO - jS plot or Brillouin diagram for the proposed circuit. This provides 
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all the information required to forecast the performance in interaction 
with electron beams, except knowledge of the actual field amplitudes 
involved, ref. Beck (1). 

The term plasma as used here denotes a partially ionized gas 
which is electrically neutral in any macroscopic volume in the absence 
of disturbances. Further, for this discussion it has been assumed the 
ions are stationary, and the electrons have no random velocities and 
suffer no collisions. The neutral gas molecules are assumed to play 
no part in the process studied. It is felt this is a reasonable model 
for the plasma. Fig. (6) shows the basic configuration of the problem. 
The helix slow wave structure has been appproximated by the sheath 
helix model. 

3. Analytical Solution. 

Trivelpiece (6), has solved the problem of a plasma filled cylindri- 
cal waveguide in an infinite axial magnetic field, starting with the 
basic Maxwell equations. He further considers the finite magnetic 
field case where he neglects the a.c. magnetic fields, (quasi-static 
approximation). He presents, for the finite magnetic field case, the 
development of a complex dielectric tensor, to take into account the 
effects of the plasma in a finite magnetic field. This is; 



12 



helix 






e 

rr 




,e re 


0 




L = € 0 


-j€ 


6r 


rt\ 

CD 

CD 


0 






0 




0 


€ 

ZZ 




C — c 




1 - 






rr 


60 = 


<o s - to 


£ 

C 



(3.1) 



(3.2) 



€ = € = ig c) ^ r 2 (3 3) 

r6 e 6r (to ) (co c 2 - to 2 ) 

. , 2 

€ = 1 — — ? (3 „ 4) 

zz to 2 

Making the quasi-static approximation, i.e. setting the a.c. magnetic 
fields to zero, we may solve the problem of matching between the helix 
and plasma column as follows: 



7 x E = -jtoB = 0 (3 o 5) 

:. e = -7 0 

This approximation is presumed good for V 2 phase < <c 2 
Now, for the equivalent dielectric, 

7*D=7" (|__ • E ) = 0 (3.6) 

which suggests the differential equation 

7 . ^ . 7 0 = 0 (3 o 7) 

The solution to which in cylindrical coordinates is: 
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In general one may assume a solution of the form 
0 = R(r) e' ,n V^ z 



(3.9) 



However, we shall restrict ourselves to the case of no variation in 0 
(i.e. n = 0). This is not essential, but will simplify the analysis. 

Also this is the case of greatest practical importance. This then leads 
to the partial differential equation 



J. 5 9^ 
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and solve 
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and now let 
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The solution to this equation, then ist 0 < r < b 
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(3.13) 



E = - A Jx (rr)e ] wt “ ^z 
r p 



T /T , , {cct-p ) 

H0 = - A — 03 ^ Jx (rr)e M z 



and for b ^ r ^ a 

E z = [bIq Or) + CKo (j3r)]e ; (t0t * ? z } 

Er = j[BIx(/3r) - CKi (jSr) j e' 1 “ ^z* (3.14) 

H0 = [BIx Or) - CKi (jSr)]e vt0t “ ^z 

and with suitable matching at the boundaries , we may reduce the above 

to; 



€ rr Jx(ra) = 

fa Jo(ra) 



(3.15) 



hifa) CKo(j8b) + FKx(jSb)]+Kx(j8a)[I 0 Q3b) - FlaQb)] 

-I 0 (/3a)[Ko(jSb) + FK 1 (j3b)]+ Ko Oa)[Io 0 b) - Fix (/5b)] 

where F is a function of the admittance Y of the sheath helix, defined 

s 

as: 

f = _J |y L (3 - 16) 

p s 

Y is the admittance of the helix. The expression for Y is given by 
s s 

Hutter (2) as: 
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Y = 

using the sheath helix approximation. It should be noted that inherent 

in use of this expression for the admittance of the sheath helix is the 

further assumption of a free space TE mode within the helix. The use 

of a TE mode in addition to the TM mode is required to match the 

boundary conditions at the helix. 

Paik, (12) has developed an analysis using the coupled mode 

approach for a helix-plasma column in the zero magnetic field case. 

The above general expression, equation (3.15) may be shown to 

reduce to that of Paik in the zero magnetic field case, (let to go to 

c 

zero) . 

Fig. 14 shows an to-/ 3 plot for the zero and infinite magnetic 

field case, with the parameters defined as follows; 

= k a cot ib 
P 

= b/a 

= k a cot ib 
c 

= pa 

= ka cot ip 
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C 3 
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y 



where k = to/C 
C = the velocity of light 
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The data was obtained by use of a digital computer. The general 

program and necessary function subroutines are given, in appendix I. 

Fig. (7) and Fig. (9) thru Fig. (17) show u; - 0 plots for various 

parameters in the general magnetic field case. Fig. (7) shows the 

principle modes, as well as a few of the higher order modes. These 

higher order modes result mathematically from the multiple roots of the 

Bessel functions, and represent higher order modes of operation. 

Since there are an infinite number of roots to the Bessel functions, 

there are an infinite set of these higher order modes in the regions 

where the propagation constant is real. A sample plot of this situation 

is shown in Fig. (8) for a typical case. The symbols used follow those 

of the computer program. Z is the right hand side of Eq. (3.15) and BB 

is the left hand side of this equation. This plot is for a fixed value of 

x. In the region where the propagation constant is imaginary, Z reduces 

to a function of two modified Bessel functions, usually represented as 

I (x) which do not exhibit this multiplicity of roots, and therefore do 
n 

not indicate higher order modes of operation in these regions. 



18 



4. Discussion. 



The cc - jS curves for the finite magnetic field case show several 
promising facets. The backward wave mode is an interesting pheno- 
menon, and indicates the possibility of a backward wave oscillator. 

Fig. (7) shows a detailed set of o.'-/3 relationships for the case 

of k a cot if) =3.0, b/a = 0.8, and k a cot if) =5.0. Fig. 8> shows 
P c 

plots of the left hand and right hand sides of Eq 3.1 5) vs ka cot if) 
for pa = 4.0. The perturbations of the helix wave as shown in Fig. (7) 
indicate that considerable coupling exists between the helix and 
plasma for this finite magnetic field case. This is encouraging, and 
indicates that one might reasonably expect to coup.e the input and 
output energy directly to and from the plasma rather than coupling to 
the electron beam and thence to the plasma as has been cone by previous 
workers, and this direct coupling to the plasma would r he r eby allow 
excitation of the backward waves as previously discussed 

The various modes of operation depicted in Fig. {"<) a r e plasma 
modes, that is, they are not space harmonic modes oi the helix, 
inasmuch as the sheath helix approximation does no: include the space 
harmonics which exist with the actual helix. Therefore in a solution 
which utilized a more exact representation for the helix one might ex- 
pect to see in addition a coupling between the plasma and higher order 
helix modes. In the limiting case as the sheath helix goes over to a 
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conducting cylinder, these modes shown would be the various radial 

modes found by Trivelpiece, (6) and denoted as the p modes where 

op 

P Q y is the Pth zero of the zero order Bessel Function of the first kind. 
From Fig. (8), we see that the solution at point a represents the per- 
turbed helix wave, the solutions to the right of point a, the modes in 
the backward wave region, and those solutions to the left of a the 
modes in the forward wave region. For example the solution at point b 
represents the first plasma wave, and this point corresponds to the 
point b on Fig. (7). The solution at point e represents the first order 
solution of the backward wave region and corresponds to the point e 
on Fig. (7). The other points correspond as they are lettered. Thus 
we may see the mathematical generation of these higher order modes, 
and may visualize them as representing higher order radial modes. There 
are of course an infinite set of these modes, only a few of which have 
been shown. 

The slight backward wave region of the perturbed heiix mode 
indicates this device might have some promise as a backward wave 
oscillator. It is noted, however, that the bandwidth of any such 
device would be relatively narrow. The forward wave characteristics 
indicate that a forward wave amplifier might be designed to be quite 
broadband, however. 

It should be noted that the higher order radial modes are probably 
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not particularly useful for interaction inasmuch as it would be difficult 
if not impossible to couple to them due to the reversal of field poten- 
tial in the radial cross section which such a mode would represent. 

Fig. (9) thru Fig. (16) show the effect of variations in oc and co 

P c 

as well as the ratio of b/a. These curves are included to provide 

some insight into the effects of a variation of these parameters. First 

it is to be noted that for to < to the lowest order plasma wave cuts 

pc 

off at to , and for to < to the plasma wave appears to cut off some- 

pep 

where between to and to . The backward waves and helix wave 
P c 

approach as a limit to when to >to , and approach as a limit to 

c c p p 



when to > to . Further, the backward waves cut off at the 
P c 






2 + to 2 



P c 

The above observations are similar to those found by Trivelpiece (6) 
in the plasma - cylinder problem. The larger b/a, the more completely 
the plasma fills the helix. It is seen that this leads to more bandwidth 
in the backward wave region. Also, the knee of the curve is less 
sharp as the helix wave goes over into a backward wave mode. This 
one might expect since a larger b/a would lead one to expect tighter 
coupling between the helix and the plasma. 

Fig. (17) shows a plot of the helix and plasma waves for the zero 
and infinite magnetic field case. The infinite magnetic field case 
shows a continuation of the effects noted previously. There appears to 
be a smooth transition to the zero magnetic field case as well. As tc 
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is decreased, the backward wave region, essentially between cc^ and 

\JtCp 3 + cc^ 2 becomes smaller and smaller, until in the limiting case 

there is no such region. The helix wave becomes more and more nearly 

like the uncoupled helix wave as to is decreased, and finally for cc =0 

c c 

we see that the helix wave shows no backward wave region, and is 

only slightly perturbed by the plasma as may be seen by comparison 

with the free helix wave also plotted on the same figure. The plasma 

wave, which at values of cc > cc appears to cut off at cc , cuts off at 

c p P 

lower and lower frequencies as cc is decreased below cc and finally 

c p 

approaches the limit as shown for the zero magnetic field or cc =0 

c 



case. 
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5. Conclusions. 



In conclusion it appears from the results of the analysis presented 
here that a plasma microwave amplifier device in which the R. F. Energy 
is coupled to the plasma directly offers some real possibilities of 
success. It appears that such a device would have several advantages 
over conventional type traveling wave devices. Principally these are 
that one might reasonable expect high power operation in a small device, 
and that one might be able to alleviate the serious mechanical and 
cooling problems inherent in higher frequency miciowave devices of the 
type presently in use. 

However, it should be emphasized that the analysis presented 
contains three basic approximations which limit the validity of the 
theory. First, the quasi - static approximation is made, which should 
be good for the low values of beam voltage used In typical tubes. 

Second, the sheath helix model is used, which is believed to be a 
reasonable model within the limitations previously mentioned. Finally, 
the assumption of a free space TE mode within the heiix is made. This 
assumption would appear to be justified since the solution obtained in 
the case of zero magnetic field agrees with Paik's analysis (12) which 
does not make this assumption. 

The dispersion relationships of Fig. (9) thru Fig. (16) do not show 
any anomalous discontinuities as were found by B. M Bulgakov et al. 
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( 13 ). 



The data contained in this paper is not sufficient in itself to draw 
firm conclusions regarding the operation of a plasma microwave ampli- 
fier device, however. The diagram does not provide any information 

regarding the magnitudes of fields involved.- and thus no insight into 
the interaction impedances involved. A theoretical study of interaction 
impedance should be an interesting and worthwhile project and is 
suggested as the subject for some future investigation. 

6. Proposed Experimental Work. 

After completing all connections for the experimental tube, it is 
proposed that the following experimental work be accomplished 

1. Make cold tests of impedance match to the helix, and deter- 
mine the losses incurred in the helix couplers over a frequency range 
of at least 1.0 to 2.5 Gc. 

2. Determine the co - j 3 relationships as indicated from the theore- 
tical work. 

4. Carry out electron beam interaction experiments. 

5. Investigate the backward wave regions shown theoretically, 
this being of special interest. These experiments should be designed 
to determine the range and magnitude of self oscillation available as 
well as bandwidth optimization by variation of field strength and 
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plasma density. 

In the experiments involving measurements of the R. F', Interaction, 
measurements should also be made to determine the impedance levels, 
gain, or loss in the helix regions, and gain or loss m the drift region. 
Further, these experiments should be carried on over a wide range of 
magnetic fields, tube temperatures, and frequencies. 

Care should be exercised to assure the cooling provided by the 
water jacket is adequate to maintain an effective temperature along 
the region so cooled. Fig. (18) shows a photograph of the laboratory 
equipment designed to work with the proposed tube. The magnet, tube, 
and heat sinks for the filament rectifiers are water cooled. Fig. 19) 
shows a block diagram of the proposed measurement setup 

The magnet is designed to provide an adjustable uniform field 
from 0 to 2000 gauss. The plasma generation filament supplies are 
designed to provide from two to six volts, and up to a maximum of 
120 amps each. The rectifiers are rated at 250 amps. , the transformer 
being the limiting element. 
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JONES GRAPH CN 8 7 

PROGR AM JONt S3 
PRINT 1 

FORMAT ( 12H TUBE DESIGN ) 

FORMAT (18H IMAGINARY GAMMA ) 
FORMAT (2X,6(2X,E15.8) ) 

FORMAT (2X,6(2X,E15.0) ) 

FORMAT ( 1 8H X , Y , HZ , HZ Y , F.-HBB ) 

FORMAT ( 5X,3( 5X,C 15.8 ) ) 

PRINT 17 

DIMENSION XI (900 > , Y1 (903 ) 
DIMENSION T ( 1000) 

Y = KACO+ PSI 
X = PA 

CONI =KP ACOT (PSI) . 

C0N2=B/ A 

F = -J(OMEGA) (EPS010N ZERO)ZP(YS) 
C0N3=5. 0 
C0N2=0 . 8 
CON 1 =3 . 0 



PRINT 7, CONI ,C0N2,C0N3 

11=50 

LL=200 

CONTINUE 

DO 40 M=LL, 15000,11 



L = M 
K A=— 1 

DO 50 K=5, 375, 1 
X=FLOATF(K )/25.0 
N = L 

AA=C0N2*X 
A0 = BE S I ( 0.0, X) 

A 1 =BES I ( 1 .0,X) 

CALL BESS (0, X,2, B0,T ) 
CALL BESS ( 1 , X,2,B1 ,T ) 
C0=BESI (0.0, AA) 

C 1 = R E S I ( 1 .0,AA) 

CALL RESS(0, AA,2,D0,T) 
CALL PE SS ( 1, AA,2,D1,T) 
DO 60 J=N, 15000, 1 
Y=FLO AT F(J)/ 2000.0 
L=J 



A=. 0000001 

I F ( SENS E SWITCH 2) 102,103 
Y=2 . *Y 1 (KK-1 )-Y 1 ( KK- 2 ) -A 
AB=Y/X 
Q= ( CON 1 / Y ) 

RA=Y **2.0 
RB=C0N3*» 2 . 0 
RC=R A-RB 

RD= ( CON 1 ) **2 . 0/RC- 1 . 0 

RR= 1 .0- ( CON 1 ) **2 .O/RC 

GAP = ABSF ( ( ( (Q)«*2.0)-1 .0)/RD) **0.5 

GAB = C0N2 *X * GAP 

XNUM1 =(AB)«« 2.0 

DENI = AO/ A 1 +B0/B1 - { B 1/ BO) *XNUM 1 
F =XNUM 1/DEN 1 

IF ( SENSE SW ITCH 1 ) 101, 120 

IF(Q-I.O) 12,13,13 
I F ( RD) 10,11,11 
IF (RD) 11,99,10 
PRINT 2 

FORMAT (16H GAMMA INFINITE ) 

CALL BES( 1*GAB,0,BESJ1»T ) 

CALL BES(0,GAB»0»8ESJ0»T ) 

• E 1 =BESJ 1 
EO =BES JO 
BB=RR«GAP*E1 /EO 
C C= C 1 » ( BO + F * ( B 1 ) ) 

DD= Dl* (AO- F * ( A 1 ) ) 
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EE=-CC*(B0+F»(B1 ) ) 
TF = DO* ( AO-r * ( A 1 ) ) 
XNUM2 = CC+DD 
DEN2 = EE +FF 
Z=XNUK2/ DEN2 
Z Y=Z/BB 

ALFA*>AC$f (ZY-1 . ) 
EPS=0. 1 



20 

10 
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IF (ALFA-EPS) 20,20,59 

PRINT 5,X,Y,Z,ZY,F,BB 

00 TO 95 

HAP = -GAP 

FO = PES I ( 0. C » GABJ 

F 1 = BESI ( 1 . 0 , GAB ) 

HCC=C1 * (RO+F *( B1 ) ) 

HDD= D1 * ( AO-F* ( A 1 ) ) 

HEE = -CO* (B0+FMB1 ) ) 
HFF = D0*( AO-F* (Al ) ) 
HEB=RR*HAP*F 1/FC 
DHN = HEE +HFF 
HNM=HCC +HDD 
HZ=HNM/DHN 
HZ Y=HZ/HPB 
H L F A = ABSF ( HZY 






• EPS=0 . 1 

IF(HLFA-EPS) 21 ,21 , 59 
59 IFISENSE SWITCH 2) 105,60 

105 IFISENSE SWITCH 3)106,107 

106 A=A+. 000005 
GO TO 102 

107 A=A-. 000005 
GO TO 102 

60 CONTINUE 

21 PRINT 15,X,Y»HZ,HZY,F,HBB 
PRINT 16 
95 KK =K -4 

XI ( KK ) = X 
Y 1 ( KK ) = Y 
NUMPT S = KK 
50 CONTINUE 

101 PRINT 5,X,Y,Z,ZY,F,RR 

PRINT S.HZ.HZY.r-.HBe.O.RI) 

CALL GRAPH ( NUMPT S , X 1 , Yl , 8 ) 

IFISENSE SWITCH 3) 12.8, 4 1 
128 11=100 
LL = M+ I I 
GO TO 43 
41 CONTINUE 
40 CONTINUE 
STOP 
END 

SUBROUTINE BESS I NO, X , KO DE , RESULT, T ) 

DIMENSION T I 1000 ) 

PI=3. 1415926536 
IF (KODE) 41,1,2 

1 CALL BES I NO, X,0, RESULT, T ) 

RETURN 

2 IF I KODE-1 ) 41,3,4 

3 CALL BES I NO , X , 1 , RESULT, T ) 

RETURN 

4 IF I KODE- 2 ) 41.5,6 

5 I F I NO- 1 ) 11,11,12 

12 PRINT 13 

13 FORMAT (34H10RDLR ERROR FOR K AND Y FUNCTIONS) 

14 STOP 14 

11 IFIX-l.) 15,15,16 
16 TI 1 ) =1 ./X 

DO 240 1=2,12 
240 TI I ) = T( I- 1 )*T( 1 ) 
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IF(NO) 21 ,20,21 

20 XGO = 1.253 3 14 1 37 3-. 156664 1816*T( 1 ) + .C801 1 1 27R2+T ( 2 )- . 09 1 3909546* 
1 T (3 ) + . 134 459 622‘3*T{4 )-. 2 2998 50 3? 0* T ( 5 )+. 3792 409 730* 1 (6) 

2-. 52 472 77331 *T ( 7 ) + . 55755 68367+T ( 8 ) - . 42 626329 1 2*T (9) 

3+ . 2 1 84 5 1 8096 *T ( 1 0 ) - . 0 668 097672 *T ( 1 1 ) + . 009 1 8933 30 *T ( 1 2 ) 

RESULT = EXPF(-t .**) /SGTTF( X)*XGO 
RETURN 

21 IF (NC-1) 12,22,12 

22 XG 1 = 1.2533141373 + . 4699 9270 1 3*T ( 1 ) - . 1 468582957+T ( 2 ) 

1 +. 1 2 80426636 *T( 3)-. 1 7364 3 1 6 37 + T ( 4 ) + . 284761 81 49*T(5 ) 

2-. 45943421 17*T(6 )+.628 338060 1 *T ( 7 ) - .66 3229 54 30 *T ( 8 ) 

3+.50SC238 5 76«T(9)-.2 5813 03 765*r< 1 0 )+. 0788000 1 1 8* T ( 11) 
4-.0108241775*T( 12) 

RESULT = EXPF(-1.*X)/SQTTF(X)*XG1 
RETURN 

15 CALL RES ( 1 , X, 1 , X I 1 , T ) 

X I 0 = T ( 1 ) 

T ( 2) =X*X/4 . 

DO 310 1=4,12,2 
310 T ( I ) =T ( 1-2 ) * T ( 2 ) 

SUM = -.57721566 +.42278420 «T(2) + .23069756 * T(4) 

1 +.03488590*T( 6) +. 00262698 *T(8) + .00010750 * T(10) 

2 + . 00000740+T ( 12) 

ANSKO = SUM - LGGF ( X / 2 . ) * XIO 

I F ( NO ) 31,30,31 

30 RESULT = ANSKO 
RETURN 

31 IF (NC-1) 12,32,12 

32 RESULT = (l./X -ANSKO + XI 1 ) / X I 0 
RETURN 

6 IF (KODE-3) 41,40,41 

41 PRINT 42 

42 FORMAT (11H1K0DE ERROR) 

43 STOP 43 

40 IF (NO-1) 50,50,12 

50 T ( 1 ) =X/3. 

DO 400 1=2,12 
400 T ( I ) =T ( I- 1 ) * T ( 1 ) 

IF ( 3 . - X ) 52,51,51 

51 SUMY = . 36746691 + .60559366 * T ( 2 ) - .74350384 « T ( 4 ) 

1 +.253C01 1 7»T(6)- .04261214 * T ( 8 ) + .00427916 * T ( 1 0 ) 

2 -.00024846*T( 12) 

CALL RES ( 1 ,X,0,XJ1 ,T) 



-.00 552740/ T (2) 
-.00072805/ T ( 5 ) 

+.00003954/ T ( 2 ) 

+ .00029333/ T ( 5 ) 



Y0 = 2./PI*L0GF( X/2. )*XJO+SLJMY 
IF (NO) 53,54,53 

54 RESULT =Y0, 

RETURN 

53 I F ( NO -1) 12,55,12 

55 R ESULT= ( X J 1 * YC-2 . / ( P I * X) )/XJO 
RETURN 

52 IF ( 4 . -X ) 59,60,60 

60 T ( 1 ) =X/3. 

DO 402 1=2,12 
402 T ( I ) =T ( 1-1 )*T( 1 ) 

SUMFO =.79788456 -.00000077 +3./X 

1 -.00009512/T(3)+. 00137237 /T(4) 

2 +.0CC14476/T(6) 

SUMS =.78539816 +.04166397 *3./X 

1 -.00262573/T(3)+. 00054125 / T ( 4 ) 

2 - . 000 1 3558 /T ( 6 ) 

YO = SUMFO/ SQRTF(X) 

IF (NO) 71,70,71 

70 RESULT =Y0 
RETURN 

71 IF (NO-1) 12,555,12 
555 CALL RES ( 1 ,X,0,XJ1 ,T) 

X JO=T( 1 ) 

GO TO 55 
59 T( 1 ) =4 . / X 



« S INF (X-SUMS) 
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DQ 401 1=2,12 

401 T ( I )=T{ I - 1 ) * T ( 1 ) 
IF (NO) 12,73, (2 
73 PO = » 3 98942279 3 



PO = .3989422793 -.0017330620 *T(2) 

-.00004876 li*T( 6 )+.30001 73565 *1(0) 
QO = -.0124669441 +.0004564324*1(2 

+ .0000 3424 68»T(6)- . 0000 1 4 2078 « T ( 0 
YO = 2./SQRTF(X)*(P0*SINF(X-PI/4. )+4./X 



•Ml)) 

* T ( 10) 



7(4 ) 
T( 10) 



72 

74 



P 1 
I 

Q 1 = 
I 

Y 1 = 



*T (4 ) 
T( 10) 
T (4 ) 
T( 10) 



+.0001754300 
-.0000037043 

-.0000869791 * 

♦ .0000442468»T ( 6 )-. 0000142078* T ( 0) +.0000032312 « 

YO = 2. /SQRTF( X) *< PO»SINF (X-PI/4 . )+4. /X*QO*C0SF ( X-PI /4. ) ) 

RESULT=YO 
RETURN 

IF (NO-1) 12,74,12 

PI = .3989422819 + . 002 92 1 8256* T ( 2 ) -.0002232030 

+.OC00580759*T(6)-.C000200920*T(8) +.0000042414 * ... 

.0374008364 - . 0006390400*T( 2 ) +.0001064741 * T(4) 

-.0000398708+T ( 61+ .0000162200*T( 8) -.0000036594 * T ( 1 0 ) 

2./SQRTF(X)«(Pl»SINF(X-3.*PI/4. ) +4 . /X*Q 1 *COS F { X- 3 . *P 1/4.1) 
RESULT=Y1 
RETURN 
END 

SUBROUTINE B ES ( NO, X, KODE , RESULT , T ) 

DIMENSION T( 100) 

KLAM=1 
KO=NO+ 1 

1009 I F ( X ) 81,101,82 
E RR= .001 

IF* ( X + ERR) 108,10 1,10 1 
IF(X-ERR) 101,101 , 108 
IF ( NO) 103, 102, 105 
T ( KQ ) = 1 . 0 
RESULT = 1 .0 
RETURN 

1F(K0) 106, 104, 105 
RESIJLT=9. 999 999999 E200 
RETURN 

105 RESULT=0.0 
RETURN 

106 PRINT 107 

107 FORMAT (55H1NEGATIVE ORDER NOT ACCEPTED IN BESSEL FUNCTION ROUTINE) 
1071 STOP 1071 

108 IF(NO) 106, 1081,1081 
1081 IF (KODE) 1 1 ,12, 1 1 

11 KLAM=KLAM+ 1 

12 J0=2*XABSF(XFIXF (X) ) 

M0=N0 



81 

82 

101 

102 



103 

104 



IF ( MO- JO ) 2,21,21 
2 MO=JC 

21 MO=MO+ 1 1 

22 T ( MO ) =0 . 

LUB=MC- 1 

T ( LUR ) = 1 .OF- 300 
GO TO ( 23,51 ) , KL AM 



23 F=2*LUR 

231 ENA ( -3 ) RAD ( MO ) 

L I L 2 ( MO ) , ENA ( T ) . 

I NA ( 1 )SAU( 233B) 

I NA ( 1 ) SAU ( 233A ) 

I NA ( 1 ) SAL ( 233A ) 

232 LDA(F)FSB(2. ) * 

STA ( F ) FDV ( X) 

233A F MU2 ( N ) FSB2 ( N ) 

233B S T A 2 ( N ) I JP2 ( 232 ) 
SUM= T ( 1 ) 

GO TO ( 39 , 69 ) , KL AM 
39 M3=M0+ 3 

DO 4 J=3»N3, 2 

4 SUM = SUM+2.*T( J) 
F= 1 . /SUM 

DO 5 J=1,K0 

5 T(J)=T(J)*F 
RESULT=T(KO) 
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69 



51 



RETURN 

K0N=7777777776777777B 

LDA(KCN),RAD(235A) 

GO TO 23 

L AC ( KCN ) > RAD (233A ) 
M3=MO+3 



UO 7 J»2 » M3 

SUM«SUM+2.*T( J) 
F= 1 . /SUM*EXPF (X) 

DO 3 J=1 , KC 



7 



8 



T(J)=T( J)*F 



RE SULT= T { KO ) 

END 

SUBROUTINE RGAM (X,Z) 

DIMENSION B ( 13) 

BID = -.422784335092 
B ( 2 ) = -.233093736365 
B( 3) = .191091 1C1 162 
B ( 4 ) = -.0245524900°7 
B ( 5 ) = -.017645242110 
B ( 6 ) = .008023278113 
B ( 7 ) = -.000804541 335 
8(8) = -.000360851491 
B ( 9 ) = .000145624324 
B( 10 ) = -.000017527917 
B ( 1 1 ) = -. 00C00262572 1 
B ( 1 2 ) = .000001323554 
B ( 1 3 ) = -.000000181220 
Z = B( 13) 

I = 12 

301 Z = Z*X + B( I ) 

I = 1-1 

IF ( 1+1 ) 302,302,301 

302 Z = (Z+1. ) *X * ( X+ 1 • ) 

END 

SUBROUTINE RCPGAMR(X.Z) 

IF (X-1.) 208,202, 203 

202 Z = 1. 

RETURN 

203 Y = 1. 

205 X = X - 1 . 

Y = Y*X 

IF (X-1.) 206,204,205 

204 Z = 1./Y 
RETURN 

206 CALL RGAM (X,Z) 

Z = Z/Y 
RETURN 

208 INT = X 
Q = INT 

IF(X-Q) 211,209,211 

209 Z = 0.0 
RETURN 

211 IF (X+1.) 212,209,213 

213 CALL RGAM(X*Z) 

RETURN 

212 Y = X 

216 X = X + 1 . 

IF (X+1.) 214,209,215 

214 Y = Y*X 
GO TO 216 

215 CALL RGAM ( X , Z ) 

Z = Y*Z 

END 

FUNCTION C E S I ( F N , Z ) 

I'F(Z) 4,5,5 
4 PRINT 6 . 

6 FORMAT! 42H BESI NOT PROGRAMMED FOR NEGATIVE ARGUMENI ) 
PRINT 7 
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7 FORMAT ( 45H BES1 COMPUTED FOR ABSOLUTE VALUE OF ARGUMENT ) 

5 Z= ABSF ( Z ) 

• I F ( MCDF ( FN , 1 .) ) 8,9,8 
9 FN=ABSF ( FN ) 

8 FK= 1 . 

FN]*FN+1 . 

PN2*PN+2. 

CALL P> C P G A MR ( F N 1 , Z A ) 

CALL RCPGAMR (FN2 ,ZB) 

A I = ( Z/2 • ) * *2 *Z B 
A = ZA +A I 
3 FK=FK+ 1 . 

A I = A I * ( Z/2. ) ** 2/ ( ( FK+FN) *FK ) 

FN2=FN2+ 1 • 

A = A + A I 

I F ( A BSF ( A I )-ABSF(A)*l0.**(-14l) 2,3,3 
2 BESI=(Z/2. )**FN* A 

END 

FUNCTION B E S K ( F N , Z ) 

FN=ABSF (FN ) 

I F ( MODF ( FN , 1 . ) ) 13,12,13 
13 F A=- F N 

B ESK = 1 . 570 79 6327* ( BESI (FA»Z)-BESI (FNjZ) ) / ( S I NF ( 3.141 592654* FN ) ) 
RETURN 
12 N=FN 

A=(-l . ) **( N+ 1 )*BESI (FN,Z )*LOGF{ 1 . 78 1 0724 1 8*Z /2 . ) 

HZ=Z/2. 

B I B A= 1 . 

B I A=0 . 

B 1 B= 0 . 

I F ( N ) 100,100, 101 

100 B=0. 

GO TO 105 

101 DO 102 KK= 1 , N 
FKK=KK 

6 IBA = B IBA/FKK 

102 B I B=B I B+ 1 . /FKK 
105. B I = B I BA 

B=BIBA*BIB 
FK= 1 . 

103 FNK=FN+FK 

B I A= B I A+ 1 . /FK 
BIB=BIB+1 ./FNK 
B I = B I* HZ **2/ ( FK* FNK) 

FK=FK+ 1 . 

TEST=BI*(BIA+BIB) 

B=B+TEST 

I F( ABSF( TEST I-AHSF (B)*10.**(-14I) 104,103,103 

104 B= ( -HZ ) **N*B 
IF(N-l) 200,201,202 

200 C=0. 

GO TO 300 

201 C=l. 

GO TO 204 

202 NA=N- 1 
NN = 1 

DO 203 J = 1 »N A 

203 NN=NN* J 
C=NN 
CB=C 

205 KD= 1 

FKD=KD 

CB= ( -1 . ) **KD*CB*HZ**2/( ( FN-FKD) *FKD) 

K D=K D+ 1 
C =C + CB 

IF(KD-N) 205,204,205 

204 C = C* HZ* * ( — N ) 

300 BESK=A+(B+C)/2. 

END 

END 
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